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We focused on elastic properties of B- and N-doped graphene in wide range of concentrations up 
to 20%. The Young's, bulk and shear moduli and Poisson's ratio have been calculated by means of 
the density functional theory for a representative set of supercells with disordered impurity patterns 
including aggregates. In contrast to earlier work, it is demonstrated that doping with nitrogen 
even strengthens the graphene layers, whereas incorporation of boron induces large structural and 
morphological changes seen in simulated STM images. Young's and shear moduli increase or decrease 
with the doping strength for nitrogen or boron, respectively, while bulk modulus and Poisson's ratio 
exhibit opposite trends. Elastic properties of samples for both types of impurities are strongly related 
to the electronic structures, especially for heavy doping (>12%). Local arrangements of dopants 
and an agregation or separation of impurities play crucial role in the determination of stiffness in 
the investigated systems. Interestingly, these findings are opossed for B- and N-contained samples. 
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I. INTRODUCTION 

Nowadays graphene, a two-dimensional monolayer 
formed out of sp 2 hybridizated carbon atoms ordered in a 
honeycomb-like lattice, attracts a lot of attention owing 
to its unique properties. Due to the hexagonal symmetry, 
its valence and conduction bands cross and have linear 
dispersion at the so called K point; these bands determine 
semimetallic character of graphene and its extremely 
high electron mobility*^ This, in connection with ex- 
cellent mechanical properties, renders graphene to be an 
ideal candidate for applications in flexible electronics 3 -^ 
and nanocomposites££ However, the ability of generat- 
ing controllable band gap in graphene is a prerequisite 
for effective applications in transistor based electronic 
devices^ Therefore, an effective functionalization that 
would open the zero energy band gap in pristine graphene 
without significant deterioration of the remaining advan- 
tageous properties is searched for, and practically, any 
thinkable way of reaching this target is currently inves- 
tigated. Substitutional doping, which refers to substitu- 
tion of carbon atoms in graphene lattice by atoms with 
different number of valence electrons such as boron or 
nitrogen, could be a recipe to achieve this goali 2 ^ 

Compared to carbon, nitrogen has one additional elec- 
tron and boron lacks one, which means that these ele- 
ments should act as electron donors and acceptors, re- 
spectively. Despite the poor environmental stability of 
chemically doped graphene^ there are many experimen- 
tal and theoretical studies that show a possibility to 
prepare p- and n-type semiconductors by substituting 
these elements into graphene JLi2.il! It has been previ- 
ously shown that, the Fermi level shifts with respect 
to the valence band top or the conduction band bot- 



tom resulting from B or N functionalization, correspond- 
ingly, and the band gap opens.— ~ n ' 15 i 16 . Boron doped 
graphene (BG) has been successfully synthesized using 
arc-discharge of graphite electrodes in the presence of H 2 , 
He and B 2 H u 16 i 17 and by the chemical vapor deposition 
(CVD) using polystyrene and HBO3.- Whereas, single- 
and a few layer nitrogen doped graphene (NG) have been 
synthesized by the CVD with NH^^ or CH 4 N 2 C£ as 
nitrogen source and using the arc-discharge method in 
the presence of H 2 and C 5 H 5 N^ or H 2 and NH 3 .±&±£ 

The BG and NG materials open various potential ap- 
plications like in graphene-based back-gate field-effect 
transistors, 9 nanosensors; 19 ' 20 fuel cells ^ lithium ion 
batterie s 22 ' 23 or even hydrogen storage^ Therefore un- 
derstanding the influence of substitutional doping, with 
boron and nitrogen, on graphene properties is proved to 
be crucial. 



In contrast to other authors, we find that the boron- 
doped graphene shows opposite elastic properties with 
respect to the nitrogen-doped graphene. These effects 
are especially well pronounced in heavily doped samples, 
not studied theoretically so far but extensively used in 
tailoring electronic devices 2 . The effects of dopants ag- 
gregation and disorder play a huge role in the determi- 
nation of Young's, bulk and shear moduli as well as the 
Poisson's ratio. In this work, we analyze many local ar- 
rangements of impurities in the graphene layer, starting 
from small concentrations and ending with 20%. There- 
fore, we extended our previous theoretical studies^ on 
electronic properties of these systems, to mechanical and 
structural behaviour at the parameters close to the tech- 
nological conditions. 



II. THEORETICAL METHODS 

Our studies of functionalized graphene layer are based 
on ab initio calculations within the framework of the spin 
polarized density functional theory (DFT)i 25 ' 26 General- 
ized gradient approximation of the exchange-correlation 
functional in Perdew-Burke-Ernzcrhof (PBE) parameter- 
ization has been applied.— Calculations have been per- 
formed using the SIESTA package. 28 ' 29 Valence electrons 
have been represented with double zeta basis sets of or- 
bitals localized on atoms, with polarization functions also 
included. The influence of core electrons has been ac- 
counted within the pseudopotential formalism. Norm- 
conserving Troullier-Martins nonlocal pseudopotentials 30 
used in our studies were cast into the Kleinman-Bylander 
separable form^ The energy cut-off for the density on 
the real space grid has been set to 800 Ry. The Bril- 
louin zone has been sampled in the 5x5x1 Monkhorst and 
Pack schemed Calculations were performed within the 
supercell scheme with the graphene layers separated by a 
distance of 80 A, large enough to eliminate any interac- 
tion. Structural optimization has been conducted using 
the conjugate gradient algorithm to achieve the residual 
forces, acting on the atoms, lower than 0.001 eV/A. 

We performed calculations for 5x5 supercells contain- 
ing 25 graphene primitive unit cells (i.e., consisting of 
50 atoms). Such supercell has been chosen in order to 
examine a wide range of dopants concentrations and a 
variety of possible distributions of the substituent atoms. 
In the described graphene supercell, one to ten B or N 
atoms have been introduced, leading to the correspond- 
ing concentration of 2-20%. In the case of two substituted 
atoms in a supercell, all eleven symmetrically nonequiva- 
lent configurations of atoms have been examined. In the 
remaining cases, twelve different randomly chosen config- 
urations have been taken into account. For each sample 
at given concentration, and with unique dopants distri- 
bution nonequivalent to other samples, the elastic prop- 
erties have been calculated and averaged at the end. The 
statistical spreads of obtained values, due to disorder and 
aggregation, are also presented. 

We simulated the scanning tunneling microscopy 
(STM) images in order to analyze the doped graphene 
structures. These images, depicted in Figure [TJ are cal- 
culated within the Tersoff-Hamann scheme^, in the con- 
stant current mode with the bias voltage (Vbias) of 0.5 
V, and visualized using WSxM code. 34 Tip shape and its 
electronic structure were not taken into account. Density 
of the tunneling current, resulting from the applied small 
bias voltage between the tip and the sample, is calculated 
according to the formula^ 
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N band 



PSTM(r,V bias ) = 

Ep-eVua, l keBZ 

where tp . ? is the wave function associated with the Ei 



eigenvalues from the energetic window set by the inte- 
grand range which contain Nf, an d states. 

We obtained the elastic properties for all considered 
concentrations and symmetrically nonequivalent config- 
urations of substituent atoms. All values for each concen- 
tration have been averaged over the investigated geomet- 
ric configurations and compared with the relevant results 
for B or N substitutions. We applied the tensile strain 
to the functionalized graphene along the one of the lat- 
tice vectors in the graphene lattice plane and calculated 
the response. The most interesting quantity, namely the 
Young's modulus, has been determined on the basis of 
the stress tensor an and the strain en components as 
follows: 



Y = — . 



(2) 



The volume of doped systems was calculated using the 
equation 



V = a-b- sin(0) • 2 



rsn s +rc(n - n s ) 



(3) 



where a and b are lattice constants, 8 is the angle be- 
tween them, re and r$ are van der Waals radii of C 
atom (equal to 0.17 nm) and substituents (equal to 0.155 
nm for N and 0.192 nm for B). The numbers n$ and n 
count dopant atoms and all the atoms in the supercell, 
respectively. The Poisson's ratio was obtained according 
to the formula 



Aa b 
~a Ab' 



(4) 



where Aa and Ab denote a change of the relevant lattice 
constant in the strained graphene layer. Bulk modulus 
(K) and shear modulus (G) (which could be easily derived 
from Eq. [5] and Eq. [3J are given, respectively, by 



K = 
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(5) 



(6) 



to complete a set of the elastic constants calculated in 
this work. 



III. RESULTS AND DISCUSSION 

A. Atomic and charge-density changes induced by 
doping 

In our previous investigations) 11 ' 12 we focused on stud- 
ies of the stability and electronic structure of B- and 
N-doped graphene. We found that these structures are 
stable, therefore it is interesting to proceed with deeper 
insight into physical and chemical properties. Important 






(a) 



ifi 



^ J*^. J* » ■*<. >*S- ^>^ _****. 








= A 



0.00 A 




(b) 



1.41 A 0.00 A 




1.09 A 



FIG. 1: (color online) Simulated STM images of graphene functionalized with B (a) and N (b) atoms are visualised above 
(Vbias=0-5V). Superimposed is a ball-and-stick model of the functionalized graphene lattice. The C atoms are depicted in 
yellow color, B in purple and N in green. 
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FIG. 2: (color online) The ball-and-stick model of four example structures of graphene functionalized boron (a,b) and nitrogen 
atoms (c,d) presented in three (5x5) supercells. In the bottom-left supercell the difference of the valence pseudocharge density 
and the superposition of spherical atomic valence pseudocharge densities (DRHO) is depicted. Concentration of substituent 
atoms: 20%. 



TABLE I: Elastic properties of pure graphene and h-BN monolayers: Young's (Y), bulk (K) and shear modulus (G) [in TPa], 
and Poisson's ratio (u), in comparison to the experimental and other theoretical values. Some papers give in-plane stiffness 
instead of Young's modulus, therefore, we have recalculated their results using as the layer thickness the Van der Waals radius 
of C atom for pure graphene or the averaged radius of N and B atoms for h-BN. 



Property 


this work 


graphene 
calc. 


exp. 


this work 


h-BN 

calc. exp. 


Y 
K 
G 

V 


1.050 
0.528 
0.449 
0.169 


0.96522 0.994 4 -£ 

0.432^ 0.445^ 0.366-0.460^ 
0.149^ 0.162a 


0.980^ l2£ 
0.582^2 

0.28C 44 


0.756 
0.443 
0.311 
0.216 


0.76922 0.781,^ 0.802ii 0.829^2 

0.461^2 
0.323^ 0.158 48 
0.2l2a 0.211*2. 0.213*2 0.217^- 



conclusion from our earlier work was that substitution 
with boron atoms is energetically more favorable than 
substitution with nitrogen atoms, which is consistent 
with other theoretical work by Berseneva et aL 36 . Wc 
have also shown how functionalization induces morpho- 
logical changes in these systems and leads to a redistribu- 
tion of the electronic charge. Moreover, incorporation of 
foreign atoms not only destroys the sp 2 hybridization of 
the carbon atoms, but also causes deformations perpen- 
dicular to the graphene sheet. Boron atoms, in compar- 
ison to nitrogen impurities, have a greater influence on 
the geometry of doped graphene layer. Mostly, it is due 
to the fact that the covalent radius of boron is larger and 
that of nitrogen is similar to the radius of carbon. We 
have found that the graphene monolayer functionalized 
with boron atoms is no longer flat, and B atoms stick 
out from the surface. Similar observation has been also 
reported by other authors , 9 i 16 i 23 and this effect is bet- 
ter pronounced for higher concentrations of substituent 
atoms. The specified changes noticed in the morphology 
of functionalized graphene, their effect on the properties, 
and finally the wide range of potential applications, en- 
couraged presented here studies of the elastic moduli. 

Before we turn to the discussion of the elastic proper- 
ties, we present shortly the structural and related elec- 
tronic properties caused by the functionalization. Scan- 
ning tunneling microscopy (STM) is a very powerful ex- 
perimental technique which allows to investigate simul- 
taneously the electronic and atomic structure of samples. 
Therefore, wc supplement our investigations with simu- 
lated STM images, which are presented in Figure [1] for 
graphene doped with boron (a) and nitrogen (b) at the 
bias voltage of 0.5 V. They provide useful information 
on the electron density of the occupied states. For high 
doping concentration of 20%, both images show deformed 
graphene structures, where impurities change the bond 
lengths and angles in lattice. We see dark holes corre- 
sponding to small current in the hexagon centers, which 
are very characteristic for the pure graphen o 18 i 37 and 
originate from destructive interference between tip and C 
atoms. Boron-doped graphene presents large bright pat- 
tern centered on the B atoms, whereas nitrogen doped 
graphene is characterized by darker areas surrounding 
substituents. In the STM images, carbon atoms are 



brighter than nitrogen because the local charge density 
of N, at the distances pictured by this technique, has 
smaller contribution to pstm than that originating from 
the first nearest neighbours of impurities. At the same 
time in B-doped graphene, the opposite situation result- 
ing from the occupied antibonding states of both boron 
and its first C-neighbours can be observed. Similar pat- 
tern was obtained by Zheng 3 ^. It is worth mentioning, 
that our simulated STM images of the N-doped graphene 
correlate well with Zhao's STM mapping. 18 

From a picture for the differential charge-density 
(DRHO) shown in Figure [5] (insets), it is clear how the 
electronic structure changes under particular arrange- 
ments of impurities. If the dopants are more clustered, 
as in Figure |2Jb) and (c), then the charge distribution is 
modified in larger areas than in the case of sparsely dis- 
tributed substituents, as in Figure [2a) and (d). Char- 
acteristic hexagons filled by negative differential charge- 
density and separated by positive differential density, 
which overlaps with bonds in most of the areas, can- 
not be clearly distinguished in the regions of clustered 
dopants. In the B-doped graphene, we observe the de- 
pletion of electron density in the B atoms vicinity and 
a transfer of the electronic charge to dopants and their 
first C-neighbors. 

Due to one additional electron in N atoms as compared 
to carbon, in the N-doped graphene the electronic trans- 
fer from dopants to the C-rings is observed. In other 
words, N substituents induce strong intervalley electron 
scattering. Charge transfer between both types of sub- 
stituents and the C atoms was also observed in some 
experiment o 16 ' 18 and theoretical works i 20 ' 35 The differ- 
ential charge-density of the doped system is much more 
locally affected by the presence of B atoms than by N 
atoms, which is consistent with the STM analysis. 



B. Elastic properties of pure graphene and 
hexagonal boron nitride (h-BN) monolayer 

As a starting point for further investigations of the 
elastic moduli of doped graphene, we present these quan- 
tities for pristine graphene and h-BN. The Young's, bulk, 
and shear moduli, and also Poisson's ratios for undoped 



systems are gathered in Table fl] 

All theoretical results obtained by us and others are 
very close. There is also a good agreement with ex- 
perimental data for Young's and bulk moduli. Discrep- 
ancy between theory and experiment is for shear moduli, 
where all theoretical values are about twice larger. 



C. Elastic properties of graphene monolayer 
functionalized with boron and nitrogen atoms 

We start with the presentation of Young's, bulk and 
shear moduli and Poisson's ratio for N- and B-doped 
graphene, in Figure [3] 

Firstly, Young's modulus decreases with the increas- 
ing number of dopants, for both types of substituents, 
as shown in Figure E^a). However, the reduction of this 
quantity in comparison to pure graphene is smaller for N- 
doped than for B-doped graphene. At the same time, the 
analysis of the geometric structure shows that boron in- 
duce more changes into the graphene layer than nitrogen 
atoms. For concentrations smaller than 14%, one can see 
that Young's modulus of NG is almost unchanged in com- 
parison to the undoped structure. Our observations are 
in good agreement with Mortazavi's works, in which sub- 
stituted nitrogen atoms do not remove flexibility of the 
graphene layer up to the concentration of 6%42, whereas 
4% of boron atoms reduce tensile strength by approxi- 
mately 8% in comparison to pure graphene 5 -^. 

Starting from the concentration of 12%, the statisti- 
cal spreads (depicted as error bars) for Young's modulus 
increase (Fig. [3] (a)), which means that the particular ar- 
rangement of dopants in graphene lattice starts to play 
an important role as far as stiffness of functionalized sys- 
tems is concerned. To visualize this, see Figure [5] which 
displays four example structures of graphene doped with 
nitrogen and boron at the highest considered concen- 
tration of 20%. For this concentration, the differences 
between Young's moduli of symmetrically nonequivalent 
configurations are equal to 0.231 TPa for NG and 0.315 
TPa for BG, respectively. Figure [H (a) and (c) refer to 
the case of the smallest Young's modulus, whereas (b) 
and (d) to the case of the highest value of this elastic 
constant observed in the investigated systems. NG de- 
picted in Figure (5Jd) is characterized by Young's mod- 
ulus higher than observed for the pure graphene. The 
stiffest structures of boron-doped graphene are those in 
which substituent atoms are aligned in the chains per- 
pendicular to the direction of the applied tensile strain. 
In case of nitrogen-doped graphene, the structures char- 
acterized by the highest value of Young's modulus are 
those in which individual nitrogen atoms are surrounded 
by carbons, including the first and the second nearest 
neighbours. In our previous investigations of stability of 
boron and nitrogen doped graphene,— we have shown 
that the configuration in which dopants are more clus- 
tered is energetically less preferable than the configura- 
tion in which dopants are separated by carbon atoms. It 



was shown by Yuge 1 ^ that creation of the B-B and N-N 
bonds is disfavored in comparison to forming of the B-C 
or N-C bonds. 

There are also differences in the charge densities be- 
tween particular arrangements of substituted atoms, 
clearly seen for higher concentrations. For a better com- 
parison, near all structures in Figure [51 we have marked 
the participation of positive values [in %] of the difference 
between the valence pseudocharge density and a super- 
position of the atomic valence-pseudocharge densities at 
presented cross-sections of the substituted systems 5 ^. In 
the case of BG, the stiffer systems have higher concentra- 
tion of positive charges in the cross-sections. Whereas in 
NG systems, Young's modulus is higher for these struc- 
tures where higher electron density is observed. 

Further, in Figure [2Kb) and (c), we present results for 
bulk and shear moduli. The magnitude of which can be 
easily calculated using Young's modulus and Poisson's 
ratio, employing formulas [5]and [6l respectively. For con- 
centration of dopants up to 16%, the bulk modulus does 
not display a significant dependency on boron or nitro- 
gen doping and its average value is almost unchanged in 
comparison to the pure graphene. In contrast for higher 
concentrations, the average bulk modulus starts to in- 
crease for BG and decrease for NG with the increasing 
number of substituents. For the dopant concentration 
of 12% and higher, the difference between symmetrically 
nonequivalent configurations starts to be pronounced. 

Shear modulus as a function of dopants, depicted in 
Fig. EKc), is very similar to the Young's modulus (Fig- 
ure [3Ka))- Shear deformations, which play an important 
role in wrinkling and rippling behavior of graphene ; 44 ! 51 
could be, like stiffness, related to the differences in lo- 
cal charge densities of doped systems. Particular ar- 
rangement of dopants in graphene layer, for concentra- 
tion higher than 10%, affects the charge carrier scatter- 
ing and the electron mobility. Shear modulus determines 
the resonance frequency of vibration modes involving tor- 
sional strains which, in contrast to the flexural strains, 
are not involved in thermoelastic loss^ It is important 
to note, that the NG-based devices are advantageous over 
the BG-based devices if one is interested in the mechan- 
ical quality factors. 

The most striking result, which distinguishes proper- 
ties of the n- and p-type doped graphene, shows up in Fig- 
ure[3Kd) which presents Poisson's ratio. This elastic prop- 
erty increases with the growing number of B atoms incor- 
porated into the graphene layer and decreases with the 
increasing number of N atoms. Graphene layers doped 
with 20% of N atoms posses average Poisson's ratio al- 
most the same as pure graphene (0.163). In contrast for 
20% of B atoms, this quantity is even higher than for 
h-BN (0.315). This ratio, which describes how easily the 
system is deformed in the direction perpendicular to the 
applied load, allows to judge the difference between boron 
and nitrogen doped graphene in this respect. For the 
same applied stretching force, the B-doped graphene will 
contract in transverse direction much more than N-doped 
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FIG. 3: Young's (a), bulk (b) and shear (c) moduli and Poisson's ratio (d) of (9,0) graphene functionalized with B (empty 
circles) and N (filled circles) as a function of the number of substituents per supercell and its percentage concentration (top 
x-axis). Symbols mean the averaged values obtained from different dopants arrangements in various supercells. The "error 
bars" are used to visualize the statistical spreads due to disorder. 
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FIG. 4: The relation between Young's modulus and Poisson's 
ratio for all the symmetrically nonequivalent configurations 
at 20% concentration of boron (a) and nitrogen (b) atoms. 
Symbols mean values obtained from disordered samples. 



graphene. For a better comparison of the differences be- 
tween both systems, we present plots of Young's modulus 



vs Poisson's ratio for the highest considered concentra- 
tion in Figure B|a) and (b). Fitting a linear function 
to each case (using the least squares method), highlights 
two opposite tendencies: the slope is negative for the 
boron-doped graphene, whereas for the nitrogen-doped 
samples, this slope is positive. 



IV. CONCLUSIONS 

In this work, we present extensive and systematic stud- 
ies of the elastic properties of B- and N-doped graphene. 
We investigated structures with 2-20% concentration of 
substituents in the graphene monolayer for various sym- 
metrically nonequivalent configurations. Our ab-initio 
studies in the framework of density functional theory in- 
clude, essential for creation of graphene-based devices, 
comparison between both type of dopants with respect 
to the concentration and particular arrangement of impu- 



rities on graphene properties. The structure and related 
electronic properties of doped graphene have been inves- 
tigated on the basis of the simulated STM images pre- 
senting samples morphology. Images of the differential 
densities between the valence pseudocharge and densities 
of the superposition of the atomic valence-pseudocharge 
(DRHO scheme) were also very helpful. We demon- 
strated that the B-doped graphene exhibits much larger 
morphological changes compared to the pure graphene 
than the N-doped samples. 

Our studies also provide valuable theoretical predic- 
tions for elastic constants (Young's, shear and bulk mod- 
uli and Poisson's ratio) of doped graphene, and evidence 
clear chemical trends, shedding light on physical mecha- 
nisms governing these effects. From the analysis of dis- 
order and aggregation, we observe a general rule that: 
aligning B atoms in the chains perpendicular to the ap- 
plied tensile strain or surrounding individual N atoms by 
C-neighbours, including first and second nearest shell, 
builds samples with the highest Young's and shear mod- 
uli. We claim that some way of incorporating N atoms in 
the graphene layer could even strengthen the doped sys- 
tem in comparison to the pure graphene, whereas all B- 
type doping noticeably reduces the stiffness. Bulk modu- 
lus and the Poisson's ratio for the concentrations greater 
than 16% and 12%, respectively, exhibit opposite trends: 
the boron atoms improve the elastic constants while the 
nitrogen atoms worsen them. In the case of N-doped 
graphene, all of the considered elastic properties remain 



almost unchanged for concentrations smaller than 14%. 

We have found that the elastic properties are strongly 
related to the electronic properties of doped systems, es- 
pecially for heavily doped samples. Particular arrange- 
ment of dopants in graphene layer and the related local 
charge-densities play important role in mechanical qual- 
ity factors in BG and NG. Each type of doping introduces 
different changes with respect to all considered proper- 
ties, and these differences cannot be neglected. 

The conclusions drawn in this work are important for 
the design of flexible electronics and nanocomposites or 
other applications in fields such as fuel cells or hydrogen 
storage, employing graphene layers. 
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